Universality in phase boundary slopes for spin glasses on self dual lattices 
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We study the effects of disorder on the slope of the disorder-temperature phase boundary near 
the Onsager point [T c = 2.269 . . .) in spin-glass models. So far, studies have focused on marginal 
or irrelevant cases of disorder. Using duality arguments, as well as exact Pfaffian techniques we 
reproduce these analytical estimates. In addition, we obtain different estimates for spin-glass models 
on hierarchical lattices where the effects of disorder are relevant. We show that the phase-boundary 
slope near the Onsager point can be used to probe for the relevance of disorder effects. 
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I. INTRODUCTION 

The critical behavior of lattice spin systems can change 
drastically when disorder is introduced. Harris [l[ 
demonstrated that whenever 2 — dv < (d the space 
dimension, v the critical exponent of the correlation 
length) the universality class of the system without dis- 
order can change when disorder is introduced Q. The 
two-dimensional Ising model Q — one of the simplest yet 
nontrivial classical spin systems — has dv — 2. When dis- 
order is added one obtains the two-dimensional Edwards- 
Anderson Ising spin glass [U, [||. However, the model is 
marginal in that the disorder only changes the singular- 
ity in the specific heat from a logarithmic to a double- 
logarithmic behavior [(| 0] • Different numerical studies 
have addressed this problem on square lattices [HHI, 
however, a detailed study of the gradual increase of the 
perturbation of the system on non-Bravais lattices re- 
mains to be performed. 

In this work we study the effects of a gradual increase of 
the disorder for small disorder perturbations for the two- 
dimensional Ising model on different lattice geometries 
and, in particular, non-Bravais lattices. We study the 
slope of the phase boundary of the disorder-temperature 
phase diagram in the limit of vanishing disorder as a 
means to detect if disorder is a relevant perturbation or 
not. 

The phase-boundary slope of the two-dimensional 
bond-diluted Ising model has been studied in detail by 
different research groups (l2l - fl4j |. In particular, Do- 
many's perturbative approach gives an "exact" result for 
the slope in the case of the bond-diluted Ising model, as 
well as the I sing model with bimodal disorder on square 
lattices [l4|, fla ]. The results stems from a simple per- 
turbation technique combined with duality arguments, 
as well as the replica method. The effective Hamilto- 
nian is separated into two parts: a nonrandom replicated 
Ising model and an operator representing the effects of 
the disorder that corresponds to interactions among dif- 



ferent replicas. The analysis, starting from this effective 
Hamiltonian, is validated when the effects of the disor- 
der are not relevant from a renormalization-group per- 
spective. Note that this approach can be applied to any 
self-dual lattice as long as one assumes that the disorder 
does not affect the nature of the critical behavior. 

Therefore, the value of the slope has certain univer- 
sal properties shared within the self-duality and the ir- 
relevant disorder. Despite nonanalyticities known to 
be present in disordered Ising models [l6|, the disorder 
seems to merely give rise to a linear shift in the criti- 
cal point (for small disorder), because most relevant op- 
erators in the disordered part of the effective Hamilto- 
nian are marginal in two dimensions. Therefore, a naive 
perturbation that evaluates only the shift of the criti- 
cal point without any further consideration can be jus- 
tified and would give an exact value of the slope of the 
phase boundary near the critical temperature in the zero- 
disorder limit. For the case of the two-dimensional Ising 
model on the square lattice, this would be the Onsager 
point [T IM = 2/ ln(l + w 2.26918 . . .] [13, E3- How- 
ever, it is unclear if this is still valid when the disorder is 
a relevant perturbation. 

In this work we compute the phase boundary slope 
for different self-dual lattices using duality arguments 
fl9l l2(il ]. We can therefore infer from the value of the 
phase boundary slope if the disorder is a relevant pertur- 
bation or not: If the obtained results do not agree with 
Domany's prediction, we can infer that the disorder is a 
relevant perturbation. Finally, we compare these results 
to the phase boundary slope of the two-dimensional Ising 
model with bimodal disorder computed using exact Pfaf- 
fian methods 12111. also validated via Monte Carlo data. 



II. HIERARCHICAL LATTICES 

Because considerably larger system sizes can be com- 
puted using hierarchical lattices, we first analyze the 



2 



phase-boundary slope on self-dual hierarchical lattices 
using the duality app roach introduced in Ref. [T^ |. Due 
to the self-duality [17| of the hierarchical lattices, the crit- 
ical point is the same as for the two-dimensional Ising 
model on the square lattice. However the relevance of 
the disorder can change depending on the geometry of 
the hierarchical lattice 123. 



A. Review: Duality analysis for spin glasses 

We first outline the duality analysis expanded to spin 
glasses. The goal is to study the phase boundary of the 
bond-diluted and bimodal Ising model in two space di- 
mensions given by the Hamiltonian 



H = -J^njSiSj. 

(ij) 



(1) 



Here Si represent Ising spins and the sum is over nearest 
neighbors. J defines the energy scale. For convenience, 
the partition function is written as a function of the cou- 
pling constant K ~ j3J (/3 = 1/T, T the temperature): 



z(k)=j2h° ktz3SiS ' 

{Si} (ij) 



(2) 



where the couplings r,j are given by 

P(Tij) = P S(r l3 - 1) + (1 - p)S(r ij + 1) (3) 

for the two-dimensional Ising model with bimodal disor- 
der and 



P(T ij )=pS(r ij -l) + (l-p)S(T ij ) 



(4) 



for the bond-diluted version of the model. For p = 1 
the standard duality analysis can be applied (l7l l23j to 
obtain the exact location of the Onsager point (without 
having to evaluate the free energy) by relating two par- 
tition functions at different temperatures, i.e., 



Z{K) = X Nb Z(K*). 



(5) 



Here Nb is the number of bonds and A = [1 + 
cxp(— 2K )]/V2 is a ratio of the local principal Boltzmann 
factors when the edge spins are parallel, i.e., xq = cxp(K) 
and its dual x^ = [exp(if) + cxp(-_ft')]/v / 2. The dual 
coupling constant is given by K* — — ln(tanhif )/2. Un- 
der the assumption that there is a single transition tem- 
perature, we identify the critical point by the condition 
K = K* = K c , i.e., the Onsager point [H H El- In 
this case the prefactor A then becomes unity. 

However, when disorder is present, a replica approach 
has to be included to determine the critical point. In this 
case the effective partition function is given by 
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where n is the replica number, and [• • -] av represents an 
average over the disorder. A similar relation as in Eq. ([5]) 
can be derived, namely 



Zn({x}) = X»°Z n ({x'}). 



(7) 



Instead of the relationship between low and high tem- 
peratures as in the case without disorder we establish 
a duality relation between the local Boltzmann factors 
Xk and x* k (§3, [25| , where k denotes the number of 
anti-parallel pairs in n-replicated bonds. The important 
local principal Boltzmann factor is the one with the edge 
spins parallel, i.e., xq = [exjp^KTij)]^; the dual one is 
xl = [(exp(KTij) + exp(-ifry-))/\/2} n ]av given by the 
2-component Fourier transform in par with standard du- 
ality procedures (23|. It follows that A„ = Xq/x . 

Because the traditional duality arguments to find a 
fixed point in an arbitrary replica number give unphysical 
results, we establish the hypothesis that the prefactor A„ 
should be unity at the critical point also for the case 
of spin glasses 0, [25| . This assumption is verified as 
correct by comparing results for different critical points 
to results from other methods. 

Recently, in Ref. (20l | an improvement of this approach 
has been proposed where one considers finite-size clusters 
to regularly sum over a subsection of spins on the square 
lattice as a real-space rcnormalization after the ordinary 
duality transformation, see Fig. [TJ This accounts better 
for the effects of the quenched randomness. Using this 
approach one obtains: 



zi s H{x^}) = (\^ 



4 S) (K (S) })> (8) 



( s) 

where N B is the number of bonds in the cluster, and 
s denotes the size of the cluster. The partition function 
is no longer written by simple two-body interactions be- 
cause of a regular summation which we denote as Z^ 
and Z*( s \ The local Boltzmann factors {x} and {x*} 
are modified by a regular summation into {x^ s '} and 
{x*( s '}, which represents these many-body interactions. 

The prefactor is also changed to An = x o^ / x< o \ which 
is the ratio of the renormalized principal Boltzmann fac- 
tors with the spins on the perimeter of the cluster fixed, 
as shown in Fig. Q] (fixed spins are represented by open 
circles). 

A single equality — X„ = 1 — then gives a more precise 
location of the critical point in comparison to the stan- 
dard duality approach because the effects of the disorder 
are properly included. 

The precision of critical points determined with the im- 
proved duality approach increases with increasing cluster 
size, similar to a renormalization-group analysis. When 
the studied lattice can be built via a recursion — as in the 
case of hierarchical lattices (2(|[27j — the real-space renor- 
malization group analysis (the regular summation over a 
part of spins in each step of the renormalization) can 
yield exact results asymptotically for the partition func- 
tion and related quantities [l9[. In this case the cluster 
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FIG. 1: Clusters used in the summation over a subset of spins 
on a square lattice for the improved duality approach. The 
full circles represent the spins summed over, while the white 
circles are treated as fixed. 
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FIG. 2: Self-dual hierarchical lattice. Full circles represent 
spins summed over while the white circles are fixed. The 
hierarchical lattice is constructed by inserting a unit cell at 
each bond. The renormalization corresponds exactly to the 
inverse process. 



size s represents the number of the step in the construc- 
tion instead of the cluster size, as shown in Fig. [2] 

(s) 

The prefactor A„ is explicitly given by the ratio of 
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[Z^{K) n } av and 
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[^ci (s) W"U where 



Z C \{K) represents the standard partition function as in 
Eq. ([2]) but limited to the cluster, and Z*(K) is its dual. 
The principal Boltzmann factors Xq 



0) 



regarded as the n-replicated partition function on the 
finite-size cluster after a configurational (disorder) aver- 
age. The quenched (replica) limit n — > of A« yields an 
equation that allows for the determination of the critical 
points of spin glasses, namely 



log z*J s \k) 



iog4 s) w 



0. 



(9) 



Starting from this equation, for instance, one can com- 
pute to high precision the location of the multicritical 
point for the bimodal (±J) Ising model (p c = 0.89082) 
using the s = 2 cluster on the square lattice. This 
also proves that there is no finitc-tcmpcrature spin-glass 
transition when the disorder distribution is symmetric 
(p = 0.5) on any self-dual lattice in two space dimen- 
sions mm. 



B. Computing phase-boundary slopes 

In the present study we estimate the value of the criti- 
cal phase boundary slope near the Onsagcr point for the 
bimodal (±J) Ising model and the bond-diluted Ising 
model. After that, we compare our results with the stan- 
dard perturbation result for the bimodal case by Domany 
for the phase boundary slope namely 14 1 



c 



1 AT 



T™ Ap 



1 



p->0 



KJa.a 



3.2091.. 



i°j/K c 



(10) 



and can be 



where {<Ji<Jj)K c is the correlator between nearest- 
neighbor spins, which is the same as the internal energy 
divided by the number of bonds. T C IM = 2/ ln(l + a/2) « 
2.26918 is the critical temperature of the two-dimensional 
Ising model without disorder [13, [U| shared within the 
self-dual hierarchical lattices. To compute the slope of 
the phase boundary close to the Onsager point we ex- 
pand Eq. © to first order for Ap <C 1 

(l - N^Ap) (logZ< s \K) log Z { S) (K)) 
+Ap£ (log Zl {s) (K) - log Z[°\K)) = 0.(11) 

Here Ng denotes the number of bonds on the finite- 
size cluster. For the bimodal (± J) Ising model, Zq and 
Zq are the partition functions without any antiferromag- 
netic interactions (or without absence of interactions for 
the bond-diluted Ising model). Z\ and Z\ are those with 
a single antifcrromagnetic interaction (or absence of in- 
teractions for the bond-diluted Ising model). The sum- 
mation (reduced configurational average for Ty ) is taken 
over all realizations with a single antiferromagnetic inter- 
action on the cluster. Next, we set the coupling constant 
as K C +AK (AK <C 1) near the Onsager point. Equation 
(fTTI) for the slope £ thus reduces to 
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FIG. 3: Unit cells of the self-dual hierarchical lattices, b is 
the length of the unit size of the hierarchical lattice. 

where we use the identity Zq(K c ) = Zq{K c ) at the On- 
sager point [see Eq. (JSJ)] . To estimate the precise value of 
the slope we need to compute the exact value of the par- 
tition function on the finite-size cluster and evaluate the 
reduced configurational average over all the possible loca- 
tions of a single antiferromagnctic interaction. The pre- 
cision of the slope computation depends on the finite-size 
cluster used. This, in turn, depends on the evaluation of 
the partition function and the reduced configuration av- 
erage, all scaling ~ O(N^) for the hierarchical lattices. 
If the calculation converges with an increasing number of 
bonds we can estimate the thermodynamic limit value to 
high precision by a simple extrapolation. 



C. Results 

On the hierarchical lattice (Fig. [3]) a simple iterative 
recursion can be used to estimate the partition func- 
tion exactly We can thus obtain an estimate 
of the slope by studying sufficiently-large clusters using 
Eq. (fT2|) . At each step of the renormalization an estimate 
is obtained for a given cluster size. Results for both bi- 
modal and bond-diluted disorder are shown in Tables U 
HI! and IIII1 respectively. The data clearly converge to a 
definite value. 

For the hierarchical lattice with 6 = 2 we sum 
up spins up to ivjj 5-17 ' 1 = 5 17 . Equation (fT2|) for 
the bimodal Ising model yields 3.209112646897937... 
for the phase boundary slope, in agreement (up 
to 13 digits) with the simple perturbative approach 
of Domany [3.209112646897908...]. For the bond- 
diluted Ising model we obtain for the phase bound- 
ary slope 1.3292579815281371 . . . using the duality ap- 
proach; the simple perturbative approach of Domany 
yields 1.3292579815281351... (agreement up to 14 dig- 
its). Summarizing, for the 6 — 2 self-dual hierarchical 
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I 

TABLE I: Results for the phase boundary slope close to the 
Onsager point for the (bimodal) ±J and bond-diluted Ising 
model on the b = 2 self-dual hierarchical lattice. The bot- 
tom line (marked with 'ED') is computed using the simple 
perturbative approach by Domany [lj, HE] . 



s 


±J 


bond diluted 


1 


3.238670311607941 


1.3313280689632047 


2 


3.213721581962695 


1.3295847287696935 


3 


3.209840157687296 


1.3293096628059718 


4 


3.209227717235686 


1.3292661587550924 


5 


3.209130853824536 


1.3292592754393175 


6 


3.209115527844001 


1.3292581862702618 


7 


3.209113102764629 


1.3292580139255745 


8 


3.209112719032133 


1.3292579866545563 


9 


3.209112658312095 


1.3292579823393163 


10 


3.209112648704037 


1.3292579816564926 


11 


3.209112647183702 


1.3292579815484458 


12 


3.209112646943131 


1.3292579815313489 


13 


3.209112646905064 


1.3292579815286436 


14 


3.209112646899041 


1.3292579815282155 


15 


3.209112646898088 


1.3292579815281478 


ED 


3.209112646897908 


1.3292579815281351 



TABLE II: Results for the phase boundary slope close to 
the Onsager point for the (bimodal) ±J and bond-diluted 
Ising model on the b — 3 self-dual hierarchical lattice. The 
bottom line (marked with 'ED') is computed using the simple 
perturbative approach by Domany ["Til , [l5( | . 



b = 3 


±J 


bond diluted 


1 


3.287620710287408 


1.3346666036511970 


2 


3.279280357814356 


1.3341234975081161 


3 


3.278702594934052 


1.3340857362686508 


1 


3.278662410902251 


1.3340831090945090 


5 


3.278659615086040 


1.3340829263029411 


6 


3.278659420560328 


1.3340829135847418 


7 


3.278659407025692 


1.3340829126998396 


8 


3.278659406083984 


1.3340829126382702 


9 


3.278659406018462 


1.3340829126339864 


10 


3.278659406013903 


1.3340829126336883 


11 


3.278659406013586 


1.3340829126336676 


ED 


3.209112646897908 


1.3292579815281351 



lattices the duality analysis yields results in perfect agree- 
ment with the simple perturbative approach. 

This implies that the disorder should be marginal or 
irrelevant on b = 2 self-dual hierarchical lattices. How- 
ever, for 6 = 3 and 6 = 4, estimates of the phase boundary 
slope using the duality approach disagree with the simple 
perturbative estimate, see Tables \H\ and ITTT1 The simple 
perturbative result is applicable to all self-dual lattices 
if we assume that the effect of the disorder operator is 
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TABLE III: Results for the phase boundary slope close to 
the Onsager point for the (bimodal) ±J and bond-diluted 
Ising model on the b — 4 self-dual hierarchical lattice. The 
bottom line (marked with 'ED') is computed using the simple 
perturbative approach by Domany [3, [H| ■ 



±J 



bond diluted 



1 

2 
3 
1 
5 
6 
7 
8 
9 

10 
11 

ED 



3.434735249924405 
3.438341429422663 
3.438917613643436 
3.439007297242741 
3.439021128918736 
3.439023254536054 
3.439023580721387 
3.439023630745810 
3.439023638415719 
3.439023639591570 
3.439023639771829 

3.209112646897908 



1.3441930333807548 
1.3444508638305390 
1.3444906297718794 
1.3444967410351193 
1.3444976788451403 
1.3444978226701897 
1.3444978447219644 
1.3444978481026600 
1.3444978486209216 
1.3444978487003697 
1.3444978487125488 

1.3292579815281351 




FIG. 4: Difference S between the estimate for the slope £ from 
the duality approach [Eq. JT2J] and the expected exact result 
[3| as a function of the cluster size Nb for 6 = 2 self-dual 
hierarchical lattices. The linear behavior in the log-log plot 
clearly illustrates the power-law convergence with the system 
size. 



irrelevant. The fact that for the 6 = 3 and 6 = 4 self-dual 
hierarchical lattices the perturbative and duality meth- 
ods disagree implies that the disorder operator of the 
effective Hamiltonian of the bimodal and bond-diluted 
Ising model is relevant. Furthermore, this suggests that 
the phase-boundary slope near the Onsager point can be 
used to probe the relevance of disorder effects. 

Finally, our results also show that the data for the 
slope £ converge for an increasing number of bonds to a 
thermodynamic limiting value as C(-^Vb) = Coo + a Ng V 
with y ps 1.15, sec Fig. @]for 6 = 2. The behavior for the 
other values of 6 is qualitatively similar. 



III. SQUARE LATTICE 
A. Numerical method 

Unfortunately, using the duality analysis on square 
lattices does not work as well as for hierarchical lat- 
tices. Finite-size effects and apparent systematic devi- 
ations prevent from an exact determination of the phase- 
boundary slope in this case. Instead of the duality analy- 
sis, we have conducted numerical simulations to compute 
a precise estimate of the slope C on the square lattice for 
p — > 1 to check the value of Domany. 

Pfaffian techniques provide access to exact partition 
functions in finite spin-glass systems, and developments 
of these techniques have become quite sophisticated, al- 
lowing exact results to be computed for larger system 
sizes and lower temperatures than are accessible with 
more traditional Monte Carlo techniques [29M3l1 |. Using 
a dual-lattice approach [3l],[32| , an algorithm has been re- 
cently developed 2l| to compute exact correlation func- 
tions of the two-dimensional Ising spin glass on the square 
lattice at all power-of-two length scales in 0(N 3 / 2 ) time 



for a system of N ■ 
C(£), is defined by 



L spins. The correlation function, 



1 



i,j s.t. nj—i 



(13) 



with r^j being the distance between sites i and j, Ng = 
Si j st r =£ 1 ^ s t ne normalization factor, and (•••) 
represents an average over the canonical ensemble. We 
can then construct a dimensionless ratio — similar to the 
Binder cumulant [33[ — for each sample 



6(4,4) 



(14) 



In the thermodynamic limit the ferromagnetic state has 
infinite correlations, i.e., the ratio tends to 6 — > 1 for 
any pair £±,£2, while in the paramagnetic state the cor- 
relations die off so that if £\ and £2 are far apart from 
one another, they will tend to 6 — > (without loss of 
generality, assuming that t\ > 1%). In order to ensure 
that £\ and £2 are far apart, it is useful for them both to 
scale with system size L. At the transition temperature 
6 approaches a constant independent of the system size 
L (up to corrections to scaling). We find empirically that 
6(L/2,L/4) has very small finite-size corrections, so we 
use this ratio for our computations. 

To extract T c from the 6 crossings, we perform a finite- 
size scaling analysis. It is expected that 



B (l 1/v (T -T c 



(15) 



where B is a unknown scaling function. The values of T c 
and v responsible for the best scaling collapse are there- 
fore our best estimates of these quantities. 
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In order to automatically extract T c and v using stan- 
dard curve-fitting techniques, we fit to the functional 
form of B. We emphasize that the form of B is unim- 
portant except as a step to automate the procedure: the 
values we extract for v and T c give a good scaling collapse 
which may be verified independently (as shown in Fig. [5]), 
and for this purpose the form of B may be discarded. We 
find that fitting to an arbitrary (four-parameter) cubic 
polynomial gives a very good scaling collapse. Two inde- 
pendent variables, L and T, appear as parameters in B. 
We therefore perform a two-dimensional curve fit using a 
nonlinear least-squares Levenberg-Marquardt algorithm. 
To correctly evaluate the statistical error bars, we per- 
form a bootstrap analysis. For each value of p and T, we 
have computed b for 10 4 different instances of disorder 
from which we built a bootstrap sample for the critical 
parameters [3J |. 

Letting q = (1 — p), simulations have been carried 
out for q = 0.0025 to 0.05 in increments of 0.025 up 
to N — I28 2 spins with 10000 disorder samples. The re- 
sults are summarized in Table IIVI and Figs. [5] and |6] In 
Fig.[5]representative scaling collapses of the raw data are 
presented. For the pure case, q = 0, there are no error 
bars because this is an exact technique. Thus the pure 
result shows the magnitude of the finite-size effects with 
this analysis. A sample result at q = 0.01 using exchange 
Monte Carlo is also presented. In this case, the Binder 
ratio 
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(16) 



is used. Here (• • • ) represents an average over Monte 
Carlo time. 31 temperatures in the range [2.0, 2.3] are 
simulated for systems up to 128 2 spins. The system is 
equilibrated for 2 16 Monte Carlo sweeps; the same time 
is used for the time average of the data. Statistical error 
bars are determined by performing a jackknifc analysis 
of 5000 independent runs. The Binder ratio is expected 
to follow the same scaling relation [Eq. (fT5")) ] as b. The 
Monte Carlo data collapse is not achieved by an addi- 
tional fit: the values of T c and v extracted for the Pfaf- 
fian data are used, illustrating that that an independent 
numerical method produces consistent results, therefore 
validating our approach using Pfaffian techniques. For 
small amounts of disorder (q < 0.01), the computed val- 
ues of v are consistent with the pure case of the Ising 
model where v = 1 [35| . and v slowly rises to a value 
of v « 1.10(4) for q = 0.05. The values of T c arc listed 
in Table IIVI Systematic errors are expected to come 
primarily from finite-size effects and therefore should be 
much smaller than the statistical errors. One can then 
extract a naive slope from each data point by using a 
finite difference with the known exact value at q = 0. 
These naive slopes are also shown in table IIVI This tech- 
nique does not give a precise estimate of C (we present a 
precise estimate from a nonlinear curve fit below) , but we 
point out that the slope is steeper than Domany's pre- 
diction for every value computed this way. This is due 




l/v 

L (T-T ) 



FIG. 5: Finite-size scaling collapse for the ratio b as a function 
of the scaling variable L 1//,y (T — T c ) for different values of q. 
Representative numerical data used for this scaling collapse 
are shown for q — 0, 0.0025, and 0.010. No fitting was done 
for q — 0: the numerical result is exact (no disorder averaging 
is necessary for q = 0), so exactly known values for T c and v 
in the pure two-dimensional Ising model are used to achieve 
the data collapse, therefore illustrating that the method used 
works. The panel on the bottom right labeled with 'MC is 
a separate finite-size scaling collapse of the Binder ratio g 
for independently-computed Monte Carlo data at q — 0.0100 
using the value of T c computed for the data in the lower left 
panel. 



to higher-order terms in the small-g expansion. 



B. Results 

Our results agree with Domany's value, confirming 
that this expansion is valid for the two-dimensional Ising 
model on a square lattice. In Fig. [S] we show the results 
of a curve fit to a third-order polynomial for the val- 
ues of T c (q). We find that linear and quadratic curve 
fits give high chi-square values and are therefore un- 
likely to accurately describe the data. We have also 
carried out curve fits for third-, fourth- and fifth-order 
polynomials: all three give similar results for the slope, 
and all three fits have similar chi-square values, imply- 
ing that the cubic fit gives correct results. The value 
of the cubic fit is the most precise of these because it 
has the fewest parameters. The obtained slope from the 
curve fit is —3. 214(17) /T c (q = 0), which is consistent 
with Domany's value, while the fit value of T c (q = 0) 
is 2.26906(16), which is consistent with the known exact 
result. 
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TABLE IV: Numerically computed values of T c (q). In ad- 
dition, a "naive" estimate of the phase boundary slope near 
the Onsager point is shown where a linear approximation is 
performed. 



q 


Tc(q) 


naive slope 


0.0025 


2.25070(11) 


3.258(19) 


0.0050 


2.23261(15) 


3.224(13) 


0.0075 


2.21349(23) 


3.273(14) 


0.0100 


2.19531(22) 


3.256(10) 


0.0125 


2.17670(59) 


3.268(21) 


0.0150 


2.15779(69) 


3.273(20) 


0.0175 


2.13834(48) 


3.295(12) 


0.0200 


2.11852(60) 


3.320(13) 


0.0225 


2.09947(42) 


3.324(8) 


0.0250 


2.08035(41) 


3.329(7) 


0.0275 


2.06058(23) 


3.343(4) 


0.0300 


2.04037(32) 


3.361(5) 


0.0325 


2.02010(35) 


3.378(5) 


0.0350 


1.99970(25) 


3.393(3) 


0.0375 


1.97927(36) 


3.407(4) 


0.0400 


1.95852(77) 


3.423(8) 


0.0425 


1.93790(34) 


3.435(3) 


0.0450 


1.91649(37) 


3.454(4) 


0.0475 


1.89517(59) 


3.470(5) 


0.0500 


1.87265(32) 


3.495(3) 



1 i 1 r 




0.01 0.02 0.03 0.04 0.05 

1 

FIG. 6: Top: Numerically-computed phase diagram for q < 
0.05. The solid line is the result of a cubic fit to the data. 
The dashed line is the q — >• limit given by Domany's value 
for the slope. Bottom: The residuals of the fit, computed by 
subtracting the value of the cubic fit from each data point, 
show no clear trends. 



IV. CONCLUDING REMARKS 

The duality analysis on hierarchical lattices shows that 
the value of the slope is not universal for all the self- 



dual lattices. This suggests that the simple perturbation 
theory for the bimodal and bond-diluted Ising model, 
without paying attention to the relevance of the disorder 
effects, can yield incorrect results. More generally, our 
results illustrate that the effects of the disorder can actu- 
ally modify the value of the slope of the paramagnetic- 
ferromagnetic phase boundary near the Onsager point 
for self-dual lattices. For the b = 2 self-dual hierarchi- 
cal lattice the value of the slope agrees with the stan- 
dard perturbative result by Domany. If the effects of the 
disorder are relevant, the value of the slope changes de- 
pending on the shape of the self-dual lattices sharing the 
same nonrandom critical point. Therefore, the value of 
the slope can be regarded as a means for classification 
of the relevance of the disorder effects on Ising models 
defined on self-dual lattices. We also numerically repro- 
duce Domany's estimate within statistical error bars for 
the square lattice using exact Pfaffian techniques. 

Several self-dual models where disorder is a relevant 
operator have a close relationship with topologically- 
protected quantum computation (36l . |37| . It is therefore 
of interest to gain a deeper understanding of the effects 
of disorder on self-dual Ising spin systems. 
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